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ABSTRACT 

Context. Relaxation theory offers a straightforward method for estimating the energy that is released when continual convective 
driving causes a magnetic field to become unstable. Thus, an upper limit to the heating caused by ensembles of nanoflaring coronal 
loops can be calculated and checked against the level of heating required to maintain observed coronal temperatures (T > 10 6 K). 
Aims. We present new results obtained from nonlinear magnetohydrodynamic (MHD) simulations of idealised coronal loops. All of 
the initial loop configurations discussed are known to be linearly kink unstable. The purpose of this work is to determine whether 
or not the simulation results agree with Taylor relaxation, which will require a modified version of relaxation theory applicable to 
unbounded field configurations. In addition, we show for two cases how the relaxation process unfolds. 

Methods. A three-dimensional (3D) MHD Lagrangian-remap code is used to simulate the evolution of a line-tied cylindrical coronal 
loop model. This model comprises three concentric layers surrounded by a potential envelope; hence, being twisted locally, each loop 
configuration is distinguished by a piecewise-constant current profile, featuring three parameters. Initially, all configurations carry 
zero-net-current fields and are in ideally unstable equilibrium. The simulation results are compared with the predictions of helicity- 
conserving relaxation theory. 

Results. For all simulations, the change in helicity is no more than 2% of the initial value; also, the numerical helicities match 
the analytically-determined values. Magnetic energy dissipation predominantly occurs via shock heating associated with magnetic 
reconnection in distributed current sheets. The energy release and final field profiles produced by the numerical simulations are in 
agreement with the predictions given by a new model of partial relaxation theory: the relaxed field is close to a linear force free state; 
however, the extent of the relaxation region is limited, while the loop undergoes some radial expansion. 

Conclusions. The results presented here support the use of partial relaxation theory, specifically, when calculating the heating-event 
distributions produced by ensembles of kink-unstable loops. The energy release increases with relaxation radius; but, once the loop 
has expanded by more than 50%, further expansion yields little more energy. We conclude that the relaxation methodology may be 
used for coronal heating studies. 
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1. Introduction 

The energy required to heat the solar corona is thought to orig- 
inate from the magnetic fields that permeate the Sun's atmo- 
sphere. The geometry of these fields is revealed by coronal 
loops, where the emitting plasma is constrained to follow the 
magnetic field: the plasma beta is extremely low (fi^O.Ol). 
Coronal loops are closed structures that emerge from the pho- 
tosphere at one location and re-enter the solar surface at another. 
The convective motions at these photospheric boundaries (or 
footpoints) are thought to reconfigure coronal fields and thereby 
cause free magnetic energy to accumulate, making the fields 
more susceptible to instability. Essentially, the kinetic energy of 
the convection zone is transported, via a Poynting flux, through 
the photosphere and stored within the coronal loop. In the case of 
a single loop, currents are created by those convective motions 
that cause the footpoints to be twisted. The approximate time 
scale for the twisting motions is long compared to the Alfven 
time, and so the coronal loop (not necessarily illuminated) tran- 
sitions through a series of force-free equilibria, which can be 
expressed as V x B = a(r)B; where r is a position vector, and 
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a = (pioj - B)/\B\ 2 is related to the parallel electric current density 
(i.e., a is a measure of the twist). 

Energy release might be triggered when a loop-like mag- 
netic field, driven by continual convective motions, reaches the 
threshold for kink instability (Hood 1992; Browning & Van der 
Linden 2003, Haynes & Arber 2007; Srivastava et al. 2010). 
Coronal magnetic fields cover many thousands of kilometers 
(L « 50 Mm) and exist within a highly conductive environment: 
therefore, in the absence of instability, magnetic fields diffuse 
slowly (fd « 10 3 yr). Nevertheless, through an ideal kink instabil- 
ity, a coronal loop may be deformed such that magnetic flux sur- 
faces are brought together. As the deformation continues, areas 
of high current are produced, allowing magnetic reconnection 
to take place. Hence, energy is released from the magnetic field 
during the nonlinear phase of an ideal kink instability, as shown 
by many three-dimensional (3D) magnetohydrodynamic (MHD) 
models (Baty & Heyvaerts 1996; Velli et al. 1997; Arber et al. 
1999; Baty 2000; Browning et al. 2008; Hood et al. 2009). 

The coronal heating problem is most pronounced within ac- 
tive regions, where the heating requirement is approximately 
10 7 erg cm -2 s" 1 , significantly higher than that necessary for the 
quiet Sun (Withbroe & Noyes 1977; Aschwanden & Acton 
2001). These regions are often observed as dense thickets of 
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transient coronal loops, the composition of which changes over 
the course of hours to days. Sudden reconfigurations are coinci- 
dent with large-scale solar flares (~ 10 34 erg). The signal strength 
of such phenomena means these events are more readily ob- 
served — detailed investigations have revealed strong evidence 
for magnetic reconnection (Fletcher 2009; Qiu 2009). Miniscule 
(nano)flares, far weaker than ones commonly observed, could 
maintain coronal temperatures if these less dramatic events oc- 
cur with sufficient frequency (Hudson 1991). This type of flaring 
might be the sort initiated by the kink instability; however, the 
amount of energy released (i.e., the difference between the en- 
ergy at instability onset and that of the relaxed field) depends 
on how much the magnetic field has altered before it relaxes. 
Bareford et al. (2010, 2011) identified the threshold for linear 
kink instability with respect to an idealised coronal loop model 
(both with and without net current). This work determined the 
subset of field configurations accessible via convective driving 
that are linearly kink unstable. One approach to calculating the 
energy release, is to represent each of these configurations within 
a nonlinear MHD code, allow the instability to take place and 
follow the reconnection dynamics until a relaxed state emerges. 
However, it is simply not feasible to run such a computationally- 
intensive process for more than a few examples. Hence, Bareford 
et al. (2010, 201 1), used relaxation theory to identify the relaxed 
states for all of the threshold (i.e., marginally unstable) configu- 
rations. 

An unstable field obeys relaxation theory if it relaxes to- 
wards the lowest energy state that conserves total magnetic ax- 
ial flux and global magnetic helicity (Taylor 1974, 1986). This 
minimum energy state is a linear force-free field; i.e., a(r) is 
invariant and Vxfi = aB. The helicity (K) indicates how in- 
tertwined a magnetic field is with itself (Berger 1999). Coronal 
loops have a non-zero normal flux at the footpoints; hence, the 
gauge-invariance of relative helicity (Berger & Field 1984; Finn 
& Antonsen 1985) makes it the more useful property: 

K = f (A + A') • (B - B ') dV, (1) 

Jv 

where A is the magnetic potential, B f is the potential field with 
the same boundary conditions and A' is the corresponding vec- 
tor potential. Helicity-conserving relaxation has been seen in 
many laboratory experiments (Heidbrink & Dang 2000; Taylor 
1986). It must be noted that helicity is subject to global resis- 
tive diffusion. Nevertheless, if localised dissipation occurs on 
small spatial scales (i.e., across shock fronts or within thin cur- 
rent sheets), the reduction in helicity will be negligible compared 
to the decrease in magnetic energy (Browning 1988, Browning 
et al. 2008). The original intention of relaxation theory was to 
explain laboratory plasma phenomena; but latterly, it has been 
frequently applied to the solar corona (Heyvaerts & Priest 1984; 
Browning et al. 1986; Vekstein et al. 1993; Zhang & Low 2003; 
Priest et al. 2005). 

The relative ease with which relaxation theory can be ap- 
plied meant that Bareford et al. (2010, 201 1) were able to gener- 
ate heating-event distributions from ensembles of idealised coro- 
nal loops, representing, albeit crudely, the population of coro- 
nal loops that exist within an active region. Each energy release 
is determined by where a loop crosses the instability threshold; 
this location is the outcome of a defined stochastic process. The 
distributions lead to an estimate for the heating rate that is just 
sufficient for coronal heating. However, the assumptions of this 
work regarding instability and relaxation theory have yet to be 
tested by a nonlinear 3D MHD code. The purpose of this paper 
is to elucidate further (Browning et al. 2008; Hood et al. 2009) 



the relaxation process and to understand how it can be applied to 
coronal loops that lack a conducting wall. We simulate a set of 
zero-net-current coronal loops that sample the linear instability 
threshold calculated by Bareford et al. (201 1). This is a larger set 
than the one investigated by Hood et al. (2009), and futhermore, 
the loops represented here feature a current-neutralisation layer 
that maintains zero net current even if the currents inside the 
loop are predominantly single signed. (Loops that carry zero net 
current are preferred since the convective motions that twist the 
loop and thereby create azimuthal field are spatially localised; 
the field outside the loop is unaffected by motions within the 
loop cross-section and therefore remains purely axial.) 

A long-standing problem has been how to apply relaxation 
theory in astrophysical contexts without the presence of con- 
ducting walls: simplistically, the relaxation should extend out to 
infinity and lead always to potential fields. Browning (1988) and 
Dixon et al. (1989) showed that relaxation theory could apply to 
volumes with free boundaries, but did not give a prediction for 
the spatial extent of the relaxed state. We find that a modification 
to Taylor relaxation theory is required before it can be used to es- 
timate the energy released by a kink-unstable loop. In contrast to 
previous work, we calculate the helicity directly at specific times 
during each simulation (Browning et al. 2008 integrated the time 
differential of the helicity in order to show the change in SK); 
thus, we are able to verify the extent of helicity conservation. 
We also examine the performance of the MHD code with regard 
to energy conservation. Any numerical dissipation will have im- 
plications for the modelling of plasma processes associated with 
heating, such as radiation. However, the plasma-/? is sufficiently 
low that any artificial resistivity should not influence the energy 
released by an unstable magnetic field, nor should it affect the 
evolution of a relaxing field. 

The paper is structured in the following manner. Section [2] 
describes the numerical code used for the nonlinear MHD simu- 
lations, along with the loop model and the equations used to cal- 
culate the magnetic field. The corresponding instability thresh- 
old for the linear kink mode is introduced, as are the positions 
of the simulated loop configurations. Section [3] presents the re- 
sults: specifically, how the different forms of energy vary over 
the course of the simulation, when the loop goes unstable, and 
how these results are affected by changes in the code parameters, 
such as spatial resolution and background resistivity. Following, 
the evolution of the loop is presented as regards magnetic field 
and current magnitude. Section [4] discusses how well the results 
fit a modified Taylor relaxation theory. Finally, in the last sec- 
tion, the results are summarised and our conclusions are given. 



2. Numerical code 

The nonlinear simulations are conducted using a 3D MHD 
Lagrangian Remap Cartesian code, called LARE3D (Arber et al. 
2001). It is written in Fortran 90 and uses the Message Passing 
Interface (MPI) to achieve parallelisation. The Lagrangian step 
uses a second-order accurate predictor-corrector step that also 
incorporates artificial viscosity, ensuring shocks are captured ac- 
curately. Van Leer (1997) gradient limiters are used at the remap 
step in order to preserve monotonicity. The divergence-free con- 
dition (V • B = 0) is maintained to machine precision by Evans & 
Hawley (1988) constrained transport. 
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LARE3D solves the resistive MHD equations given by 

-V-(pv), (2) 



dp 

dt 



d 1 

— (pv) = -V-(pvv) + —(VxB)xB - VP + Vo-, (3) 

— = Vx (vxB) - Vx 77 , 

dt v ' \ fx ) 



dt 



-V-(pev) - PV-v + 7// + ecr, 



with specific energy density, 



6 = 



(r-i)p 



(4) 



(5) 



(6) 



where p is the mass density, v is the plasma velocity, B the mag- 
netic field, P the thermal pressure, 77 is the resistivity (not mag- 
netic diffusivity), / is the current density, y - 5/3 is the ratio of 
specific heats, and jiq is the magnetic permeability. Viscous heat- 
ing is represented by the last term of equation ([5]), which also in- 
corporates artificial viscosity (Wilkins 1980) in order to capture 
the heating effect of shocks. This heating term is expressed as 
the product of the rate of strain tensor, 



Udvi 
2\dj 



dvj_ 
di 



and the shock tensor, 



= p/( 



v\c f + v 2 



l\s\) 



(7) 



(8) 



where Cf is the fast magnetoacoustic speed, / is the distance 
across a grid cell in the direction normal to the shock front, s 
is a similarly localised strain rate (the subscripts i and j denote 
the different spatial coordinates), and the artificial viscosity co- 
efficients vi = 0.1 and V2 = 0.5 are constants. The form of equa- 
tion ([5]) is derived from the Rankine-Hugoniot jump conditions 
and the values of the coefficients have been chosen such that nu- 
merical oscillations behind shock fronts are prevented. Note, the 
force equation ^ also acquires a viscous term. Gravitational ef- 
fects are ignored in this study, as are thermal conduction and ra- 
diation. The simulations are concerned with how the magnetic 
field changes in response to the kink instability; specifically, 
how much magnetic energy is released and how the field subse- 
quently evolves. Conduction becomes important some time after 
the energy release and later, radiation is the dominant process. 
Note, numerical studies have shown that conduction can act on 
MHD time scales (Botha et al. 2011): the amount of energy re- 
leased from the field is unaffected, but the kinetic energy parallel 
to the field is much reduced. 

The MHD equations are made dimensionless by replacing 
the variables with dimensionless equivalents. For example, 

r* p* B* 

r = — , p = — , B = — , 
Po B\ 

where asterisks denote dimensional variables, R& is the loop ra- 
dius, po the initial mass density, and B\ the initial axial field at 
r = 0. The other variables are expressed in a similar manner; 



L 



R~, 



f_ 



V 

va 



p* 



where £* = 20 R& is the loop length, t^ = R& / va is the ra- 
dial Alfven transit time, va=B\/ y/ioPo the Alfven speed, and 
Po = BI/jjq the magnetic pressure. The specific energy density, 
current density and resistivity (6, / and 77) also have reference 
variables that can be expressed in terms of po and B\ : 



^0 



B\_ 
j^oPo 



Jo 



Bi 



7/0 = ^RsVa- 



Values appropriate for a coronal active region can be ob- 
tained by setting R b = 1 Mm, p = 1.6726 x 10" 13 kgm" 3 and 
Bi = 50 G. Hence, the length becomes 20 Mm, v A « 10 Mm s" 1 
and 7/ ~47TX 10 6 Dm. 

The resistivity is taken to be non-uniform in these simula- 
tions, 



n = 

V = rib + , 



\J\< J evil , 
m>/crit, 



where 77b is the background resistivity (normally set to zero, 
since actual coronal resistivities are approximately yUom 2 s -1 ) 
and 7/ c = 0.001 is the anomalous resistivity, which is only 



switched on when the current reaches or exceeds J C1 



15. The 



value of / cr it is set so that it is significantly higher than the maxi- 
mum current at the start of the simulation. Super-critical currents 
appear when, during the nonlinear evolution of the kink instabil- 
ity, current sheets begin to form and decrease in thickness. The 
anomalous resistivity is intended to capture the dissipation oc- 
curring at scales below the grid resolution: at this scale, resistiv- 
ity is enhanced by small-scale plasma instabilities. 

The computational domain is a 3D staggered grid: physical 
variables are not calculated at the same place for each cell in the 
domain. This approach improves numerical stability and allows 
conservation laws to be included in the computation. The do- 
main size is L x = Ly = 6 (-3: +3) and L z = 20 (-10: +10). Initially, 
the loop axis follows the z-axis and the loop radius is r = 1 ; there- 
fore, the simulated loops all have an aspect ratio of 20. The loop 
is line-tied at z = - 10, +10, which means, at those z-coordinates, 
the velocity components are set to zero. The velocity compo- 
nents are zero at the boundaries for all directions. The normal 
derivatives of magnetic field, energy and density are zero at all 
boundaries. The simulations are run with two grid resolutions: 
128 2 x 256 (low) and 256 2 x 512 (high). It is assumed that a re- 
sult is not a numerical artefact if it is consistent across both res- 
olutions. 



2.1. Initial configuration 

Some previous studies have used LARE3D to simulate the ap- 
plication of kink perturbations to a straightened line-tied coro- 
nal loop, see Gerrard et al. (2002), Gerrard & Hood (2003), 
Browning et al. (2008) and Hood et al. (2009). The initial equi- 
librium model used in the latter two studies was extended by 
Bareford et al. (2011) to include an outer current-neutralising 
layer so as to ensure the loop has (at least initially) zero net 
current: this improves on the model used by Browning & Van 
der Linden (2003) and Bareford et al. (2010), which allowed 
loops to have net current (i.e., an azimuthal field was usually 
present in the potential envelope). All currents are now created 
by convective motions local to the loop footpoints. Hence, a cur- 
rent neutralisation layer is introduced here, defined such that the 
azimuthal field (B e ) always falls to zero at the loop boundary 
(Re)\ therefore, B e is zero in the potential envelope. The loop's 
radial Of-profile is approximated by a piecewise-constant func- 
tion featuring three parameters (Figure [T]): the ratio of current 
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photosphere 



o 



cross section 




photosphere 



Fig. 1. Schematic of a straightened coronal loop in the x-z plane (left) 
and in the x-y plane (right). The loop, comprises a core (dark grey), 
an outer layer (light grey) and a current neutralisation layer (blue); the 
whole loop is embedded in a rectangular potential envelope. The core 
radius is half the loop radius {Ri'.R^.Re'.R'B = 0.5:0.9:1:3, where Ry is 
the distance from the initial loop axis to the edge of the envelope). The 
loop aspect ratio (L/R 6 ) in this figure is 20. 



layer current) is found, for a given (a\, a 2 ), by numerical solu- 
tion of B 3 q(R&) = 0, ensuring that the net current is zero and that 
the azimuthal field vanishes outside the loop, see Eq. ( [T4| ). From 
the nondimensionalisation of the magnetic field, the field coeffi- 
cient at the core is B\ = 1. 

The linear kink instability threshold for this current- 
neutralised loop was determined by the CILTS code (Van der 
Linden 1991; Browning & Van der Linden 2003) — it uses a 
bicubic Hermite finite element method to calculate the growth 
rates and eigenfunctions for specific line-tied ^-configurations. 
In contrast to the stability space for a loop of net current 



4 
2 

a 2 
-2 
-4 
-6 



unstable 



stable (mixed B ) 



stable 



stable (mixed B v ) 



unstable 



to magnetic field is a\ in the core, a 2 in the outer layer, a 3 in 
the neutralisation layer and zero in the potential envelope. The 
free parameters are oc\ and a 2 , whereas a 3 is dependent on the 
first two and is determined by the requirement of zero net cur- 
rent. The magnetic field is continuous everywhere, whereas the 
current has discontinuities, and the outer surface of the poten- 
tial envelope, representing the background corona, is placed at 
R<b = 3 (three times the loop radius); this is far enough away that 
the boundary conditions do not influence the plasma evolution. 

The fields are expressed in terms of the well-known Bessel 
function model, generalised to the concentric layer geome- 
try (Melrose et al. 1994; Browning & Van der Linden 2003; 
Browning et al. 2008). The field equations for the four regions 
(core, outer layer, neutralisation layer and envelope) are as fol- 
lows: 



(9) 

< r < R x , (10) 
(11) 

Ri<r<R 2 , (12) 



B lz = BiJoQadr), 
B w = a^Ma^r), 

B 2z = B 2 J (\a 2 \r) + C 2 Y (\a 2 \r) , 
B 2e = (r 2 (B 2 Ji(\a 2 \r) + C 2 Y 1 (\a 2 \r)), 

B 3z = B 3 J (\a 3 \r) + C 3 Y (\a 3 \r) , (13) 
B 3e = o- 3 (B 3 Ji(\a3\r) + C 3 Fi(|a 3 |r)) , R 2 <r<R b , (14) 

Ba z = B 4 , (15) 
B 4e = 0, R& < r < R$ , (16) 

where <x; = ^ (i = 1, 2, 3) represents the sign of o^-. The fields 
must be continuous at the inner radial boundaries, R\,R 2 and R&. 
(The positions are R\ = 0.5, R 2 = 0.9 and R& = 1, so that most of 
the loop is similar to the one described by Bareford et al. (2010), 
but with a thin current neutralisation layer between R 2 and R&.) 
Therefore, the coefficients Bj and Cj (j = 2, 3, 4) are determined 
by the requirement of continuity of the magnetic field at all inter- 
faces (Bareford et al. 2011). The value of a 3 (the neutralisation 



«1 

Fig. 2. The linear kink instability thresholds for LjR b = 20. These 
thresholds have been cropped by a pair of dashed lines that indicate 
where B z (r) starts to acquire a mixed polarity. The right threshold is 
sampled by a selection of five marginally unstable configurations (black 
circles). 



(Bareford et al. 2010), the instability threshold is open, see 
Figure [2] This is very similar to the threshold found for a close- 
fitting conducting shell (Browning & Van der Linden 2003), 
since in the case of zero net current, perturbations quickly fall 
to zero beyond the loop radius. A consequence of the Bessel 
function model is that the axial field changes sign if the values 
of a\ and a 2 are too large. This reversal is unphysical and hence, 
introduces a restriction to the stable region of Figure |2j namely, 
all stable configurations should have B z (r) of uniform sign. The 
new stability space is closed by excluding those field profiles that 
have a B z (r) of mixed polarity, since this could not be achieved 
directly by footpoint motions of an initially unidirectional field. 
The filled circles of Figure [2]identify the loop configurations (see 
also Table [T} that will be simulated by the LARE3D code (the 
initial field profiles for Loops B and D are given in Figure[3]). All 
of these configurations are unstable to the ideal kink instability. 

Loops of uniform twist (loops A-C) are a more likely result 
of correlated convective driving (Bareford et al. 2011), where 
the threshold is approached via a series of steps with da\ « 6a 2 
— this is the reason why the part of the threshold curve where 
a 2 > is more finely sampled compared to a 2 < 0. Note, the sec- 
tion of threshold on the left of Figure |2]is merely the negative of 
the section on the right; hence, it is not necessary to sample both 
threshold sections. 

The configurations listed in Table [T] encompass two types of 
loop; one class where oc\ and a 2 have the same sign (A-C) and 
the other where these parameters have opposing signs (D and E). 
Loop B has been chosen as representative of the first loop type 
and Loop D of the second. Henceforth, Loop B will be referred 
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Table 1. The a-profiles, axial fluxes and field coefficients for the simulated loops. 



Loop 


ari 


a 2 


a 3 




B 2 


c 2 


B 3 


c 3 


B 4 


A 


2.42 


2.4 


-13.08 


15.3 


0.1 


-0.0079 


2.42 


-0.59 


0.55 


B (Ami) 


2.25 


1.5 


-8.71 


18.1 


0.79 


-0.2 


-0.27 


2.37 


0.64 


C 


2.15 


0.53 


-4.95 


20.8 


0.61 


-0.15 


-0.94 


-1.85 


0.74 


D (L mix ) 


2.54 


-1.0 


-0.84 


21.0 


0.91 


0.50 


0.92 


0.38 


0.74 


E 


2.8 


-2.7 


3.82 


20.4 


0.26 


1.32 


-1.79 


0.026 


0.72 


Stable 


0.5 


0.1 


-0.96 


27.8 


0.97 


-0.0078 


1.21 


0.63 


0.98 



1.2 
1 

0.8 
0.6 
0.4 
0.2 


-0.2 
-0.4 
-0.6 



Loop L u 



Loop L n 



Stable Loop 



-1 o 




Fig. 3. The initial analytically-determined B z (solid) and B e (dashed) profiles for 
Table[TJ are also shown. 



and L mix . The field profiles for the stable loop (last row of 



to as L un i (uniform twist) and Loop D will be denoted by L m [ x 
(mixed twist). 

Each of the loops indicated in Figure [2] is subjected to a dis- 
turbance of the form, 



c\e 



-4r 4 



cos) — z)cos (kz - 6) 



(17) 



where v r is the radial component of the perturbed velocity, L 
is the loop length, r = ^x 2 + y 2 is the radial coordinate, £=1.1 
is the wave number, 6= arctan(x/y), and the constant c\ = 0.01 
reduces the amplitude so that the perturbation is initially linear. 
This disturbance initiates the kink instability. 

Both L un i and L mix are simulated for low and high resolu- 
tions; only the high resolution is used for the other loop configu- 
rations. Each simulation runs until the magnetic field appears to 
have settled into a lower energy state. 

3. Numerical results 

3.1. Energy and resistivity 

Loop L un i (a i = 2.25, - 1.5) is linearly kink unstable, and the 
numerical simulation (Figure [4j middle row) shows that it is also 
nonlinearly unstable. The magnetic energy, W, has an initial (di- 
mensionless) value of 155.3 when integrated over the entire sim- 
ulation domain, 



W 



the internal (U) and kinetic (E) energies are also volume inte- 
grals, 



U 



y 



^-jJpdV, E=^Jpv 2 dV. 



The initial magnetic energy integrated over the loop L un i only is 
much less (Ws ~ 20). All magnetic energies can be dimension- 
alised by making a simple correction to equation (27) given in 



Bareford et al. (2010): since here the axial flux is not normalised 
to 1, the dimensionalising multiplier must first be divided by if/^, 
which is the square of the non-dimensional axial flux over the ra- 
dius R<b (Table[T] column 5). For loop L un i, this gives a multiplier 
of 4.8 x 10 19 (assuming a radius of 1 Mm and a mean axial field 
strength of 50 G): thus, the dimensionalised initial loop energy 
is roughly 10 21 J. 

At the onset of instability, W undergoes a decrease, coinci- 
dent with a rise in U and with a much more modest increase 
in kinetic energy (£ max ~ 0.2). The maximum current, / max , also 
rises just before the decrease in W, indicating the formation of 
a helical current sheet. The nonlinear instability starts at around 
t = 50 tA and within the next 50 1^, approximately 70% of the to- 
tal energy release has been achieved. Magnetic energy reduces 
more slowly after t- 100 ^a- The release of magnetic energy is 
of a similar size for both resolutions, and significantly, larger 
currents are recorded at high resolution, which is expected; oth- 
erwise, spatially-confined changes in current are missed and 
anomalous resistivity is reduced. The fact that the maximum cur- 
rent increases with resolution is indicative of current sheet for- 
mation. These structures have (possibly) infinite current density, 
so higher resolutions should reveal larger values of the maxi- 
mum. 

In the top two cases of Figure |4j ideal MHD and zero back- 
ground resisitivity, it can be seen that there is no detectable 
change in energy until around 40 tp, — this is the time when 
a current sheet starts to form. The adjective ideal is italicised 
because the rate of magnetic diffusion in the corona is easily ex- 
ceeded by numerical resistivity; therefore, any reduction in mag- 
netic energy associated with a stable configuration is artificial. 
Hence, for the stable case, in which no current sheets ever form, 
there is no sudden onset of dissipation, only gradual dissipation 
as expected. In the bottom row of Figure |4](77 b = 10" 4 ), there is a 
continual dissipation of magnetic energy due to Ohmic effects; 
however, even here, there is a clear onset of enhanced dissipation 
at the point of current sheet formation. This all agrees with pre- 
vious work (Browning et al. 2008; Hood et al. 2009) and clearly 
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Fig. 4. Loop L uni : the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current for ideal 
MHD (top row), for resistive MHD with 77b = (middle row), and for 77b = 10~ 4 (bottom row). The initial magnetic energy has been subtracted from 
the magnetic energy plots (solid lines, left column). The critical current (/ crit = 15) is indicated by the grey dashed horizontal lines (right column). 
For the 77b = case, the maximum current plots are from the high (black) and low (grey) resolution simulations. 



indicates that energy dissipation (conversion from magnetic en- 
ergy to internal energy) is closely associated with the existence 
of current sheets. 

However, Figure [4] also shows that energy is not conserved 
for the two cases mentioned (ideal MHD and 77b = 0). At high 
resolution, only ~ 68% of the energy released from the mag- 
netic field is converted to internal energy, and by the end of the 
simulation, less than one percent of the energy release is in the 
form of kinetic energy. Halving the spatial resolution worsens 
the energy conservation by almost 10%. This loss of energy is 
due to spatially unresolved current sheets, resulting in numerical 
diffusion. The results for Ideal MHD (i.e., zero magnetic dif- 
fusivity) are very similar to the ones produced for the resistive 
MHD case with no background resistivity. This shows that, for 
these cases, Ohmic dissipation does not contribute significantly 
to the rise in internal energy, which is instead mainly caused by 
viscous heating. The use of a non-zero background resistivity 
(/7b = 10" 4 ) changes the plots. The reduction in magnetic energy 
starts right away and is more drawn out and the maximum cur- 
rent is much less than when 77b = 0. These differences are all con- 
sistent with an increased resistivity. Crucially, Ohmic heating is 
now clearly present and, with the viscous heating, accounts for 
nearly all of the dissipated magnetic energy. In terms of energy 
conservation, LARE3D performs better when 77b > 0: the inter- 
nal energy increase is now approximately 96% of the magnetic 
energy decrease (the final kinetic energy is less than 0.2% of 
6W). A background resistivity of 10" 4 appears to be the small- 
est value that effectively minimises the effect of the numerical 
resistivity, such that the results are likely to be a reasonable de- 



scription of the thermodynamics. If 77b is halved, the percentage 
of the magnetic energy release lost to the simulation increases to 
around 10%. Further tests have shown that energy conservation 
is not improved by a lowering the critical current and thereby 
causing the anomalous resistivity (t7 c ) to be applied earlier in the 
simulation. 

A sufficiently large background resistivity does substantially 
mitigate the amount of energy lost by numerical resistivity, at 
least for L un i; but 77b = 10" 4 is not realistic by coronal standards, 
and Figure]?] (bottom left) reveals a significant drop in field en- 
ergy before the kink instability takes effect. This initial decline 
is caused by global Ohmic diffusion rather than magnetic recon- 
nection, since current sheets have not yet formed. The numerical 
resistivity comes about when current sheet widths begin to fall 
below the grid resolution; although magnetic energy continues to 
be dissipated, it is not converted to other forms of energy (e.g., 
internal or kinetic). However, shocks can still be resolved during 
the unstable phase, and these contribute to the internal energy 
via viscous heating. If the background resistivity is large enough 
it will limit currents and thereby prevent current sheets from be- 
coming too thin. Numerical (that is to say artificial) resistivity 
will be a factor during the conversion process when 775 = 0: ap- 
proximately 32% of the reduction in magnetic field energy is not 
accounted for by the final internal energy — this is also true for 
the other loop configurations. However, virtually all of this arti- 
ficial resistivity is occurring during the energy conversion. The 
drop in magnetic energy is a robust result. To demonstrate fur- 
ther, we have used ideal MHD to simulate a loop (a\ = 0.5 and 
(X2 = 0.1, see Table [T] bottom row) that is well within the stable 
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Fig. 5. Loop L mix : the temporal variation in energy (magnetic, internal and kinetic), heating (Ohmic and viscous) and maximum current (low and 
high resolution) for resistive MHD with r/ h = 0. The different plot lines follow the same scheme as that used for Figure [4] 



region (as shown in Figure [2]). In the absence of an instability 
(and any applied resistivity), the energy declines by around one 
thousandth of one percent over 300 tp^\ this reduction is three or- 
ders of magnitude smaller than the energy release caused by the 
kink instability. 

Loop L m i x (a i = 2.54, - -1.0) has a core that is oppositely 
twisted with respect to its outer layer. Figure [5] shows the same 
correspondence between the magnetic and internal energies that 
was seen for the previous loop. Again, the magnetic energy re- 
lease is of the same size for both resolutions and, at the higher 
resolution, significantly larger currents are recorded; although 
the size of the release is around half that found for L um . Note, 
the initial magnetic energy is higher than the value given for the 
other loop. This is a consequence of setting B\ = l: L m i x is less 
twisted and therefore has a higher axial flux, which results in a 
lower dimensionalised magnetic energy. 

The general trends for magnetic, internal and kinetic ener- 
gies are consistent between resolutions (for both loops), and 
most importantly, so is the size of the magnetic energy re- 
lease. Therefore, simulations at higher resolution are not re- 
quired — this paper will proceed with results taken only from 
the 256 2 x512 simulations. The results for the other loop con- 
figurations on the threshold (Figure [2]) also suggest that linear 
instability evolves to a nonlinear stage, which gives rise to cur- 
rent sheets, magnetic reconnection and most importantly, shock 
heating. The following sections will present results based on 
a current-dependent resistivity and zero background resistivity 
— this is more compatible with the coronal environment. The 
breakdown of energy conservation associated with this parame- 
ter choice can be ignored since we are only interested in how the 
magnetic field behaves during and after the instability. However, 
this issue will have to be addressed for detailed studies of the 
thermodynamic evolution. 

3.2. Magnetic field and current magnitude 

Now we examine the magnetic field (and critical current dis- 
tribution) at specific times during the simulations. Figure [6] 
shows how field lines, originally located within the core, be- 
come kinked as the instability takes hold. The dark grey field 
lines are drawn from the bottom boundary (or footpoint) and the 
light grey field lines are drawn from similar locations at the up- 
per boundary. Loops L un i and L m i x follow the same course of 
events. Initially, the field lines are intertwined; then, during the 
growth of the instability, the currents become critical (indicated 
by the yellow, orange and red areas) and anomalous resitivity 
is applied. The dissipative effects of this anomalous resistivity 
have a minimal contribution to the internal energy (Section [3T| ); 
instead magnetic energy is dissipated by the application of vis- 



cous heating at shock fronts. Hence, we also show a proxy for 
viscous heating, | <x* |, which is equivalent to the shock tensor 
divided by p L (vi Cf + V2L\s\) — \cr*\ is represented by the 
cyan, blue and purple colours. In general, shocks are coincident 
with current sheets: i.e., the areas of high viscous heating, as 
indicated by | cr* |, are cospatial with the largest currents. Note, 
the kink instability is more pronounced for L un i than for Lmix. In 
terms of azimuthal field, L m i x is the weaker of the two (Figure 
[3J, however, some importance should be attached to the fact that 
is twisted both ways. Opposing twists appear to mitigate 
the growth of the instability and thereby limit the energy release. 
If we increase the values of a\ and but keep the opposite 
signs, such that the total azimuthal field strength is comparable 
to L un i (e.g., Loop E), the energy release increases only slightly 

Figure [7] shows cross sections of L un i at z = (halfway along 
the loop) and L m i x at z - - 2, which is roughly the centre of 
the only patch of significant viscous heating for Lmix at t - 60 
see bottom row of Figure [6] Again, the colours represent current 
magnitudes (left column) and viscous heating (right column), 
and the plot times are the same as those used for Figure |6j i.e., 
shortly after the start of the kink instability. At this time, cur- 
rent sheets of narrow width start to form; furthermore, nearly all 
of these current sheets are associated with shock formation and 
viscous heating. 

The unstable phase is over quickly (At « 50 t&) and by the 
end of the simulation the (reconnected) field lines have straight- 
ened considerably, indicative of a low constant-a configuration. 
The areas where shocks have formed must be dispersed through- 
out the loop volume in order for helicity to be more evenly re- 
distributed and thereby create a linear Qf-profile. The reduction in 
the azimuthal components of the field lines should cause a radial 
expansion of the loop. At the initial equilibrium, the inward ten- 
sion force of the azimuthal field is balanced by the outward mag- 
netic pressure of the axial field; thus, if the tension decreases, the 
loop must expand before equilibrium can be regained. This be- 
haviour is clearly demonstrated by Figure [8j note the change in 
scale for the x and y axes. The expansion of the loop is mostly 
associated with the reconnection of initially twisted field lines 
inside the loop with the ambient axial field. In Figure [5j the fi- 
nal state calculated by the numerical simulation is overlaid with 
magnetic field vectors in the x-y plane, which are consistent with 
a cylindrical configuration, bounded by a current-neutralising 
layer: the arrows follow each other and the arrow sizes initially 
increase away from the axis, and then diminish before the loop 
edge. 

Loop Lmix expands less than L un i, which is possibly due to 
the fact that the latter releases more energy. A single weakly- 
twisted flux tube best describes the final state for both loops. The 
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Fig. 6. Magnetic field lines originating from the bottom left footpoint (dark grey) and from the upper right footpoint (light grey) are shown at 
t = 60 t A for L uni (top) and L mix (bottom). At the onset of instability, two plots are shown: one with isosurfaces of current (left) and the other with 
isosurfaces of | <r* |, a proxy for viscous heating (right). 
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Fig. 7. The spatial variation of current magnitude (left) and a viscous heating proxy (right) across the loop cross section at the apex (i.e., where 
z = 0) for L uni (top) and at z = - 2 for L mix (bottom). 
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Fig. 8. The spatial variation of current magnitude across the loop cross section at the final time of the simulation for L uni (left) and L mix (right). 
All currents are now well below the critical value. The arrows are magnetic field vectors. 



following sections will show that the properties of these loops 
are consistent with relaxation theory. 



3.3. Helicity conservation 

DeVore (2000) showed how to calculate the magnetic helicity 
over an entire coronal volume above a photospheric bounding 
surface. The first step is to work out the magnetic vector poten- 
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tial for a current- free field that has the same distribution of verti- 
cal magnetic flux at the lower boundary. DeVore begins by deriv- 
ing an expression for the scalar potential, using Green's function 
for Laplace's equation as the integration kernel, 



Table 2. Helicity at three times during the simulations (rj^ = 0) 
for all five kink-unstable loops. 



c(x,y,z?,t) 



+3 , B z (x',y, -10, t) 

dy' zV , — , (18) 

3 r' 



where r' - ^/(x - x f ) 2 + (y - y') 2 + (z') 2 - The grid domain used 
by LARE3D has a Cartesian geometry: the coronal loop is ini- 
tially represented as a straight cylinder within a rectangular box. 
The x and y axes extend between -3 and +3; hence, the integral 
limits given above. The photospheric boundaries are located at 
the limits of the z axis (z = -10, +10) and z = is the loop apex; 
Eq. ( [T8] ) uses the first boundary position. The vector potential is 
constructed using, 

X+io 
dz (f> c (x,y,z',t); (19) 

which becomes, 

A c (x,y,z,t) = — I dz I dx' I dy' 
2tt J z J_3 J_ 3 

l[( x -^)S-(y-y)i] (20) 



7z vy-3 J-3 

B z (x\y, -io, Of 



•M3 



(r') : 

when the derivatives are moved inside the integral. Now the 
gauge-invariant vector potential can be specified as, 

A(x,y,z,t) = A c (x,y,z,t) - zx f dz B(x,y,z',t), (21) 

J-io 

by subtracting the helicity due to the potential field, and Eq. pT) 
can be re-expressed by expanding the cross product of the sec- 
ond term, 

A(x,y,z,t) = A c (x,y,z,t) + I dz 

J-io 

x [B y (x,y,z\i)x - B x {x,y,z\t)y\. (22) 
Finally, the gauge-invariant magnetic helicity is 

K = J A BdV 

X+3 yo+3 
dx I rfy I A(x,y,z,t)-B(x,y,z,t). (23) 

3 J-3 J-10 

The geometry used by DeVore differs significantly from that 
used here (Figure [T}, which features two separate photospheric 
boundaries at the limits of the z axis. Fortunately, the relative po- 
sitions of the two boundaries mean that if the flux is cancelled 
for one it will be cancelled for the other, and so the lower bound 
z coordinate can simply be set to -10. 

Equations ([T8])-([23]) have been implemented, using the five- 
point Newton-Cotes integration formula. The marginally unsta- 
ble loops (Figure [2]) have zero net current initially and should 
continue to do so during the simulation; thus, outside the loop 
the helicity is zero. This means the helicity, calculated using a 
straightforward cylindrical geometry, can be compared easily to 
that calculated for a Cartesian geometry, where the loop is en- 
closed within a rectangular box. The helicity is zero everywhere 
in the additional volume between the surface of the rectangular 
box and the outer edge of a cylindrical potential envelope. 



Loop 


Initial 


Instability 


Final 


AK/AW 


A 


12.3 


12.26 (f = 50f A ) 


12.22 


0.09 


B (Ami) 


12.29 


12.27 (f = 60f A ) 


12.28 


0.03 


C 


10.47 


10.46 (f=100f A ) 


10.5 


0.19 


D (Anix) 


6.12 


6.12(f = 60f A ) 


6.11 


0.13 


E 


1.16 


1.14(f = 50f A ) 


1.18 


1.32 



Table [2] gives the helicities for each loop at three different times. 
The second time is the time of instability; i.e., when the loop is 
furthest from equilibrium. The helicity appears to be conserved: 
it varies little over the course of the simulation. For loops A-D, 
the helicity varies by less than one percent; whereas for Loop E 
the helicity increases by ~ 2%. 

Next, we calculate the ratio of the helicity variation to the 
change in magnetic energy, 



AK 
AW 



ln(£final/ initial) 



ln(W fina l/Winitial) 



(24) 



both AK and AW are weighted by initial value. For four of the 
five loops, AK is around one order of magnitude lower than 
AW (Table [2] fifth column); these results are comparable with 
Browning et al. (2008). The relationship AK <<c AW implies that 
magnetic energy dissipation is taking place on small spatial 
scales (i.e., shock fronts). One of the loops (E) has AK>AW; 
however, it is no coincidence that this loop also has the lowest 
helicity (initial = 1.16). The coarseness of the grid sampling used 
to calculate Eq. (23 ) — the x and y dimensions are sampled at 
every fourth cell and the z dimension at every other cell — is 
too great for such a small helicity. Hence, for this last loop, the 
grid sampling is increased from 4x4x2 to 2x2x2, which gives 
AK/AW « 0.44. A finer sampling and/or higher resolution is re- 
quired to bring this ratio down to below 0.1. 



4. Partial relaxation model 

4.1. Analytical calculation 

In Browning & Van der Linden (2003), Browning et al. (2008) 
and Bareford et al. (2010), a loop with net current was assumed 
to relax such that it expanded to fill the entire potential enve- 
lope: the Qf-profile was invariant between r = and r = 3 (R®). 
The relaxed alpha was identified by assuming that if/ (axial flux) 
and K/i// 2 (the normalised helicity) were conserved over the loop 
and envelope, in accordance with Taylor relaxation. The only 
limit to the relaxation is the position of an (unphysical) bound- 
ing wall. Hence, the relaxed state always represented a threefold 
radial expansion of the initial state. Later, Bareford et al. (201 1), 
considered a range of relaxation radii: Re<r<R^. Some form 
of partial relaxation is more likely to be relevant to the zero-net- 
current case; it is known from simulation results, presented here 
and in Hood et al. (2009), that reconnection is of limited extent, 
leaving much of the external field undisturbed. This contrasts 
with previous results for loops carrying net current (Browning et 
al. 2008), in which the disturbances in the nonlinear phase of the 
kink generally extended to the boundary of the simulation. 

Bareford et al. (2011) maintained zero net current by fix- 
ing a neutralising loop surface at a specified radius (r = R[). The 
neutralising surface is imposed by fixing the field coefficients of 
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Fig. 9. The unstable and relaxed loop. K/if/ 2 is conserved over the re- 
gion enclosed by the dashed circle, which is the outer boundary (R ( ) of 
the relaxed loop. 



where Wrfan , an) is the energy of the threshold state and Wc(ac) 
is the energy of the relaxed state. 



Relaxation Alpha 



Energy Release 




Fig. 10. The variation with relaxation radius (R[) of relaxed alpha (or/-, 
left) and of dimensionless energy release (SW, right) for L uni (solid line) 
and for L mix (dashed line). 



the potential envelope, so that they do not change during relax- 
ation. In the relaxed state, the envelope is the region between Rt 
and Ry. Naively, it might have been expected that the partially- 
relaxed state would consist of a constant a field embedded di- 
rectly in a potential field (a = 0) with no layer of reversed cur- 
rent. However, such an arrangement does not match the simula- 
tion results, see Section 4.1.1 Furthermore, it would be unphysi- 



cal for the partially -relaxed state to include a region of azimuthal 
field extending to infinity, which is the consequence of allowing 
net current, since in the initial state, the azimthal field is zero 
outside the loop. The existence of current sheets (here broadened 
into a narrow current-neutralising layer) bounding localised re- 
laxed states has also been noted by Gimblett et al. (2006). The 
axial flux is conserved such that if/t of the threshold state is equal 
to 07 of the relaxed state. The subscript C denotes the relaxation 
radius, Rt', it is the radial upper bound over which the associ- 
ated property is calculated — the lower bound being the axis. 
For example, if/t is the axial flux from r = 0to r = R( (similarly, 
KrB is the helicity from r = to r = R$, the outer edge of the po- 
tential envelope). This method of conservation is illustrated by 
Figure [9] Left, is the marginally unstable threshold state; the ra- 
dius that this loop will expand to is indicated by the dashed cir- 
cle; right, is the expanded relaxed loop. We might expect that, ei- 
ther the relaxed loop flux (if/t) matches the initial loop flux (if/ 6), 
or it matches the initial flux within the radius (Rt) eventually at- 
tained. The former is correct if the field freely expands into the 
surrounding field. The latter choice is made if the loop radially 
expands by reconnection; i.e., the loop eats into the surrounding 
axial field. This outcome a grees m uch better with the simulation 
results (again see Section 4.1.1| ). If the loop did not reconnect 
with its surroundings and expand to a radius of 1.8 7^^ (as for 
Luni), then the axial field would drop by a factor of 1.8 2 , which is 
not observed. (The slight mismatch between the numerical and 
analytical axial field profiles seen in Figure [TT] is probably due 
to some limited free expansion of the loop.) 

The relaxation alpha (ac) is determined by ensuring that the 
axial flux and helicity integrated over the dashed circle in Figure 
[9] (left), match the values obtained when these same proper- 
ties are integrated over 0<r<R[ of the relaxed loop (Figure [9] 
right). Also, helicity is absent from the potential envelope sur- 
rounding a zero-net-current loop; hence, K& = K$ (unlike mag- 
netic energy, W& ± W$). Thus, we do not need to choose which 
value of helicity to use. Once at is known, the energy release 
can be determined analytically; 



The relaxation radius (Rt) remains as a free parameter, al- 
though sensible values can be inferred from the numerical re- 
sults. Figure 10 shows how at and SW vary with relaxation ra- 
dius. As expected, these figures show an inverse relationship be- 
tween at and SW. In general, SW increases with Rt, however, 
this relationship is not linear; beyond a moderate expansion of 
50% (Rt = 1.5) the energy release is -99% of its maximum for 
L m i x and -80% for L un i. This property of diminishing returns 
is also true for the other three loops indicated in Figure [2] This 
means that the energy release is insensitive to the choice of re- 
laxation radius, so long as it is assumed that Rt > 1.5. 



4.1.1. Numerical Relaxed States 

The final numerically-determined state is merely the last snap- 
shot provided by the simulation; it is expected to be close to the 
analytically-determined relaxed state. Each loop is simulated for 
at least 300 tp,\ so, the sooner the instability occurs, the closer the 
loop will be to a fully relaxed state by the end of the simulation. 
The Cartesian components of the analytical relaxed field are as 
follows, 



B x (r) = -B e (r)(y/r), 
B y (r) = B e (r)(x/r), 
B z (r) = B.JoQatlr), 



(26) 
(27) 
(28) 



SW = Wt(a iu a i2 )-Wt(at), 



(25) 



where B e (r) = (ad\at\) B\J\(\at\r) and r<R[. We generate re- 
laxed configurations for every value of Rt between 0.9 and 3.0 in 
increments of 0.01. Then, for each field component, it is checked 
which of the 211 possibilities has the lowest chi-squared value 
when compared with the numerical values for the same com- 
ponent. These field comparisons (B x , B y and B z ) are performed 
over the x dimension for a selection of fifteen y-z coordinate 
pairs (y <e {-1, -0.5, 0, +0.5, +1}, z e {-5, 0, +5}) — the horizon- 
tal dashed lines of Figure [13] (left) indicate the y positions of the 
field plots. 

As a result of the kink instability, the loop axis will be shifted 
from the origin in the x-y plane. The new axis position will be z- 
dependent and can be recalculated by assuming that the relaxed 
loop is current neutralised (Figure [5]): moving inwards from the 
edge of the envelope, the loop boundary is detected whenever the 
a- value rises above some minimal value. Any axis shift is taken 
into account before the analytical results are compared with the 
numerical data. Figure JJ_ presents a subset of the analytical- 
numerical comparisons for L un i at y = 0. 

In general, the analytical plots, such as the ones in Figure 
1 1 show a good agreement with the numerics; however, the Rt 
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Fig. 11. Loop L uni : a comparison between the B x (top row), B y (middle row) and B z (bottom row) magnetic field profiles obtained numerically 
(red line) and analytically (black line). The latter is calculated from the a c and R c that best fit the numerical plot, which is taken from the final 
frame (t = 400 t A ) of the high resolution LARE3D simulation (rj h = 0). The comparisons are done at y = for different z coordinates, z = - 5 (left 
column), (middle column) and 5 (right column). 



and at that best fit the red lines are independently chosen for 
each of the forty-five subplots (there are fifteen y-z positions and 
three field components). Figure 12 (left) shows the variation in 
these best-fit parameters: each symbol, by virtue of its position, 
gives the Rf-ai pair that best fits the numerical data extracted at 
specific y and z coordinates within the simulation domain — the 
symbol type denotes which field component is being matched. 
The use of two x-axes (one for R(, the other ac) means that all 
forty-five symbols can be clearly distinguished (the y-axis has no 
meaning beyond a simple sorting of the data points). The derived 
averages are (R C ) = 1.16 ± 0.29 and (a c ) = 0.2% ± 0.31; these 
two properties have a one-to-one mapping. The dimensionless 
energy release is 5.87 ± 1.11. Despite the large scatter for ac, 
the deviation for SW is comparatively modest, this is because 
d(6W)/d(R[) is small when Rc « 1.8 (Figure [TO) right). Figure 
12 suggests that the final numerically-calculated state is tending 
towards a relaxed loop that can be characterised by a localised 
invariant a-profile. (Note, all points would lie on a single ver- 
tical line if there were a perfect match between the analytical 
and numerical models.) The outlying points on the left have ac 
values that are further from {ac) than those on the right, since 
the relationship between ac and Rc is not linear (Figure [T0| left) 
- this explains the large deviation for the relaxation alpha. Low 
levels of magnetic field that fluctuate around the zero line tend to 
be fitted by low ac values; whereas values significantly greater 
than the mean imply that high currents still exist within the fi- 
nal numerical state. Both of these effects can be ameliorated 
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Fig. 12. Loop L uni : the best-fit R[ - a t pairs when r] h = (left) and when 
z/b = 10~ 4 (right). Black circles are for those best fits determined from 
B z profiles, red plus signs are for B x and blue crosses B y . 



by re-running the simulation with background resistivity; once 



again, 7/ b = 10 . Figure 



12 



(right) yields (R c ) = 1.74 ± 0.11 and 
(ac) = 0.23 ± 0.09. The mean dimensionless energy release be- 
comes (SW) = 6.06 ± 0.32. The resistivity smooths out low-level 
noise and restricts the current, and thereby reduces the deviation 
associated with the analytical fit. 
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Table 3. Analytical and numerical comparison for the kink- 
unstable loops A-E. 



j-B/(\j\\B\) 



Analytical 



Numerical 



The overall impressive level of agreement demonstrated for 
Luni, also extends to other positions along the instability thresh- 
old, see Table [3] This table also confirms that the analytically- 
calculated helicities of Loops A-E are in approximate agree- 
ment with the numerical values, which were derived accord- 



Simulation 


Rc 


at 


K 


\SW\ 


K 


\SW\ 


A 


1.9 


0.25 


12.89 


8.47 + 0.42 


12.19 


8.31 


B (Ami) 


1.76 


0.28 


12.92 


5.87 + 1.11 


12.3 


5.8 


C 


1.8 


0.2 


11.03 


3.7 + 0.42 


10.52 


3.15 


D (Anix) 


1.59 


0.19 


6.45 


2.38 + 0.28 


6.14 


2.15 


E 


1.18 


0.19 


1.22 


2.63 + 0.11 


1.15 


2.58 



ing to the procedure discussed in Section |3.3| The correspon- 
dence between the numerical and analytical energy releases is 
the most significant finding. There is evidence to suggest that 
this correlation persists even when different settings are used for 
the LARE3D parameters controlling resistive MHD (Section[3T 
and Figure]?]). In addition, these results are consistent with pre- 
vious work (Browning et al. 2008; Hood et al. 2009). 

4.2. Final Current Distribution 

The jaggedness of the numerical field profiles shown in Figure 
[TT| clearly indicates that there are deviations from a simple 
locally-constant a profile. We show the value of a computed 
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j3 = 2P/\B\ 2 
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Fig. 13. The value of a (left) and plasma-p (right) over the x-y plane 
for Ami (?7b = 0). The grey circle approximates the loop cross section at 
z = 0and£ = 400? A - 

over the midplane (z = 0) of L un i, see Figure [13] A grey cir- 
cle based on the current magnitude plot of Figure [8] (left) has 
been used to approximate the shape of the cross section. There 
are many, albeit confined, areas that are far from the calculated 
mean, at = 0.28. The influence of the initial (the a- value for 
the current neutralisation layer) can be seen in the limits of the 
colour bar. In addition, plasma- J3 (Figure 13 right) has, com- 
pared to initial values of 10" 3 , undergone localised increases of 
up to two orders of magnitude. 

The cosine of the angle between the current and the mag- 
netic field is plotted in Figure 14 note, positions outside the 
loop are left unplotted. The current is either parallel (white) 
or anti-parallel (black) to the field for a high percentage of 
the cross sectional area (positions that are far from parallel, 
20° < jlB < 160°, account for approximately 30% of the loop 
cross section). Almost three quarters of the energy released from 
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Fig. 14. Left, the cosine of the angle between the current and the mag- 
netic field over the x-y plane. Right, the cosine of the angle between 
the magnetic field and the gradient of the pressure over the same area. 
Positions outside the grey circle, representing the loop cross section, 
have not been assigned a colour. The background resistivity is zero. 

the field becomes internal energy (i.e., thermal pressure), which 
may be associated with departures from a force-free state (al- 
though Taylor relaxation would predict a uniform pressure dis- 
tribution). However, it seems that the loop is heading towards 
an approximate balance of forces (jxB = VP): Figure 14 (right) 
shows that the angle between the field and the pressure gradient 
is on average close to zero. Interestingly, the parallel and anti- 
parallel areas are often found next to each other and therefore 
might be expected to cancel at later simulation times. The plots 
of Figures [13] and [14] are qualitatively similar to those taken at 
other z coordinates that are not too close to the footpoints (e.g., 
-5<z<5). 



5. Summary and conclusions 

A nonlinear 3D MHD code has been used to simulate the evo- 
lution of a set of zero-net-current cylindrical loops. These loop 
configurations have been identified by a linear analysis as be- 
ing marginally kink unstable (Bareford et al. 2011). The sim- 
ulations show that the instability quickly enters a nonlinear 
phase and magnetic energy declines sharply before leveling off. 
Furthermore, the amount of energy released matches the amount 
predicted by Taylor relaxation (Table [3]), taking account of the 
fact that the relaxation is localised. Evidence for helicity conser- 
vation was presented, and the change in helicity was shown to be 
much smaller than the drop in magnetic energy. The implication 
of this result is that energy diffusion is occuring on small scales 
compared to the global length scale; i.e., within shocks associ- 
ated with magnetic reconnection. The low values of kinetic en- 
ergy during the unstable phase (compared to the internal energy) 
imply that these shocks occur near to reconnection sites. (We 
find that dissipation within current sheets only becomes signifi- 
cant when rjb = 10" 4 .) Further research could reveal exactly how 
magnetic reconnection spreads through the loop volume in re- 
sponse to a kink instability, which would also reveal where the 
loop is heated and when. It is this widespread dispersal of recon- 
nection and shock heating that ensures helicity conservation. At 
present, we cannot confirm the nature of the shocks: slow-mode 
shocks are expected since only shocks of this type can reduce 
magnetic field strength. 

Relaxation theory also predicts that the final relaxed state 
should have a constant ^-profile. Although the final numerical a- 
profile still retains much fine structure, the final magnetic fields 
are well-modelled by a (localised) constant-a profile with some 
fluctuations superposed — this suggests the fine-scale structure 
is self-cancelling (i.e., it integrates out). Furthermore, the en- 
ergy of the final numerical state is very well matched by the 
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energy of the same constant-a state (the field energy is insen- 
sitive to the spikeness of the numerical data). The property of 
zero net current is retained after the instability. Typically, the 
loop expands radially, the field reconnecting with that present in 
the potential envelope. These results justify the choices made 
by Bareford et al. (2011) regarding the details of the relax- 
ation process. Nevertheless, we were concerned that the thin- 
ness of the current neutralisation layer (Figure [T]) might have 
influenced the results. Hence, we also ran a resistive MHD sim- 
ulation (with r/b = 0) for a zero-net-current equilibrium that pos- 
sessed a smooth ^-profile, see Case 3 of Hood et al. (2009). The 
constant parameter, A= 1.62, was set such that the equilibrium 
was on the threshold of instability. We found that the simula- 
tion results were similar to those produced by L m i x . Numerically, 
the energy release was 1.2, which was again consistent with the 
analytically-determined value. 

It appears that the assumption of a Taylor-relaxed state, sub- 
sequent to a kink instability, has been verified by the work pre- 
sented here. The relaxation does not extend over the full numer- 
ical volume, but over a region of smaller extent (out to a ra- 
dius R[, which is less than the full radius, R<b). In this sense, 
the relaxation is partial. A relaxed state can only be identified if 
the relaxation radius is known; at present, it is unclear how R[ 
can be precisely determined from the field configuration at insta- 
bility onset. However, the analytical work has revealed that for 
marginally-unstable loops, the energy release varies little with 
relaxation radius once R[>\.5 (Figure 10); hence, a calculation 
of the energy release does not necessarily require a precise pre- 
diction for R[. 

Energy release could be limited if the unstable loop attains 
an equilibrium that is less than fully relaxed (i.e., the ^-profile 
remains nonlinear) and still conserves helicity. There is perhaps, 
for some field configurations, another constraint that decides the 
relaxed state, such as the topological degree of the field line map- 
ping between the ends of the loop, as investigated by Yeates et 
al. (2010). They examined two braided magnetic field config- 
urations (one based on the simple pigtail braid and the other 
more complex). Both configurations underwent turbulent relax- 
ation, leading to a final state that conserved topological degree 
and was less relaxed than that predicted by Taylor theory — the 
final state for the pigtail braid featured two flux tubes of oppo- 
site twist. Nevertheless, it is possible for the Taylor-relaxed state 
and the state that preserves topological degree to coincide. In our 
case the invariants given by Yeates et al. do not provide any extra 
constraint, making our results consistent with their predictions. 
This would explain the level of agreement between the LARE3D 
simulations and Taylor relaxation. 

An issue for further research concerns the interaction of con- 
vective driving with the relaxation process. The LARE3D code 
could be used to help resolve this issue. It should be possible 
to choose a loop configuration (i.e., a set of a-parameters) that 
is just inside the threshold for linear kink instability and then, 
trigger the instability by applying a pre-determined velocity 
profile (vq) at one of the footpoints. Loop curvature has also 
not been considered. The linear stability codes require an 
analytical form for the magnetic fields. If a loop is to retain its 
curvature, it can only be simulated numerically, which means 
choices have to be made concerning loop parameters (e.g., 
length, radius and a-profile). Usefully, those straightened loop 
configurations that are kink unstable and are likely to be reached 
by convective driving have been identified. These configurations 
could be adapted to include curvature and re-simulated within 
LARE3D. This would reveal what effect, if any, curvature has 
on the energy release precipitated by kink instability. Of course, 
this procedure could also be applied to other improvements; 



e.g., gravity (with p(z)), conduction and radiation. However, a 
feature that improves the realism of the loop model may not be 
important as regards kink instability and Taylor relaxation. 
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